The first time you use a new version of R you will have to install all your packages. Try installing the raster and rgdal packages if you don’t have them already:
## Loading required package: raster
## Loading required package: sp
## Loading required package: rgdal
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Loading required package: lwgeom
## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: osmdata
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
Now it is installed in your computer but isn’t loaded You will have to ‘import’ your library package every time you start up R.
library(raster)
library(rgdal)
Raster data is stored in a format similar to that of a matrix, except you have added data about georeferencing.
Let’s read in one raster and take a look
notice that the working directory doesn’t need to be a path directly below your files of interest:
CWD_apr = raster('./data/raster/CWD/cwd2010apr.tif')
Get info about the raster object
CWD_apr
## class : RasterLayer
## dimensions : 1120, 872, 976640 (nrow, ncol, ncell)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## source : cwd2010apr.tif
## names : cwd2010apr
## values : 0, 150.9694 (min, max)
Plot the image:
plot(CWD_apr)
there are a variety of different plot options. Here are a few important ones
plot(CWD_apr,main = 'Title')
plot(CWD_apr,col= heat.colors(10), main='Water Deficit' ) # 10 designates the number of hex color codes that get generated
Now use ?plot to change col to a different preset color scheme
?plot
## Help on topic 'plot' was found in the following packages:
##
## Package Library
## terra /cloud/lib/x86_64-pc-linux-gnu-library/4.2
## raster /cloud/lib/x86_64-pc-linux-gnu-library/4.2
## graphics /opt/R/4.2.2/lib/R/library
## base /opt/R/4.2.2/lib/R/library
## sf /cloud/lib/x86_64-pc-linux-gnu-library/4.2
##
##
## Using the first match ...
Rasters like matrices can make basic math functions very easy. Note that rasters can be used in basic operations like CWD_apr *2 or CWD_apr^2.
However teo rasters with different extents, cell sizes, or projections CAN’T be used. You would need to reprojected and aligned first.
CWD_apr_sq = CWD_apr^2
plot(CWD_apr_sq, main='CWD Squared')
Now lets add two different rasters. Use the raster() function to read in the CWD for May
#CWD_may = raster('fill/in/the/path/file.tif')
CWD_A_M = CWD_apr + CWD_may
plot(CWD_A_M)
We can also directly set the cell values of rasters. For instance, it
is often useful to have a raster filled with a fixed value, such as
0. In this example we will use the setValues()
function to update the values of CWD_may to all be 0, and
store it in a variable called zeros.
zeros = setValues(CWD_may,0)
plot(zeros)
Note that zeros has the same extent, resolution, and
projection as CWD_may
print(zeros)
## class : RasterLayer
## dimensions : 1120, 872, 976640 (nrow, ncol, ncell)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## source : memory
## names : cwd2010may
## values : 0, 0 (min, max)
print(CWD_may)
## class : RasterLayer
## dimensions : 1120, 872, 976640 (nrow, ncol, ncell)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## source : cwd2010may.tif
## names : cwd2010may
We can also use boolean queries to update certain values. Here we
update certain values of CWD_may to 5000 where
CWD_may>150 is TRUE with the following
syntax:
CWD_may[CWD_may>150] = 5000
plot(CWD_may)
Note that we have now changed the values of
CWD_may… if we
want to original values again, we will need to read them back in again
with
CWD_may = raster('./data/raster/CWD/cwd2010may.tif').
In this section we will try to get the average water deficit for
California for the whole year of 2010. To do this we will need to
isolate the files we are interested in using - those ending with
.tif.
First lets get a list of all files in our raster folder:
list.files(path="./data/raster/CWD/")
## [1] "cwd2010apr.tfw" "cwd2010apr.tif" "cwd2010apr.tif.aux.xml"
## [4] "cwd2010apr.tif.xml" "cwd2010aug.tfw" "cwd2010aug.tif"
## [7] "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml" "cwd2010dec.tfw"
## [10] "cwd2010dec.tif" "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"
## [13] "cwd2010feb.tfw" "cwd2010feb.tif" "cwd2010feb.tif.aux.xml"
## [16] "cwd2010feb.tif.xml" "cwd2010jan.tfw" "cwd2010jan.tif"
## [19] "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml" "cwd2010jul.tfw"
## [22] "cwd2010jul.tif" "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"
## [25] "cwd2010jun.tfw" "cwd2010jun.tif" "cwd2010jun.tif.aux.xml"
## [28] "cwd2010jun.tif.xml" "cwd2010mar.tfw" "cwd2010mar.tif"
## [31] "cwd2010may.tif" "cwd2010nov.tif" "cwd2010oct.tif"
## [34] "cwd2010sep.tif"
Notice that .png and twf xml files are also listed. We can narrow this down by using a grep style pattern search.
Here we will only list files that have .tif in the
name:
list.files(path="./data/raster/CWD/",pattern = '.tif')
## [1] "cwd2010apr.tif" "cwd2010apr.tif.aux.xml" "cwd2010apr.tif.xml"
## [4] "cwd2010aug.tif" "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml"
## [7] "cwd2010dec.tif" "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"
## [10] "cwd2010feb.tif" "cwd2010feb.tif.aux.xml" "cwd2010feb.tif.xml"
## [13] "cwd2010jan.tif" "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml"
## [16] "cwd2010jul.tif" "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"
## [19] "cwd2010jun.tif" "cwd2010jun.tif.aux.xml" "cwd2010jun.tif.xml"
## [22] "cwd2010mar.tif" "cwd2010may.tif" "cwd2010nov.tif"
## [25] "cwd2010oct.tif" "cwd2010sep.tif"
almost there… we still have problems with .tif.* files like tif.xml
or tif.aux.xml. We can limit these using what I will call an
ANTI-wildcard. “$”. '.tif\$' will tell the computer the
file needs to END in .tif. Let’s also store the path/file
name by setting full.names=TRUE
list.files(path="./data/raster/CWD/",pattern = '.tif$', full.names=TRUE)
## [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
## [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
## [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
## [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
## [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"
Let’s store that information:
CWDS = list.files(path="./data/raster/CWD/", pattern = 'tif$', full.names=TRUE)
Challenge: In teams read in and sum all of these rasters, then divide by the total number of files.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
Finally, let’s write out our results:
# writeRaster(x=yourraster,filename="./output/CWD_mn_2010.tif")
R also has the capactiy to handle multi-band or multi-layer raster
‘stacks’. In this
case we will be using a stack to combine multiple files of the same type
(CWD) over time (months of 2010). As you will see there are multiple
advantages to handling data in stacks.
The only requirement is that all the images must have the same: - Image extent - Resolution - Projection
In this case our CWD data is already identical in shape and
projection. All we need to do is use a list of files like
CWDS that we created earlier.
CWDS
## [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
## [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
## [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
## [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
## [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"
and create a stack from the data:
cwd_stack = stack(CWDS)
Now looking at the properties we can see that cwd_stack
has 12 layers.
cwd_stack
## class : RasterStack
## dimensions : 1120, 872, 976640, 12 (nrow, ncol, ncell, nlayers)
## resolution : 1080, 1080 (x, y)
## extent : -375034, 566726, -616284.9, 593315.1 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
## names : cwd2010apr, cwd2010aug, cwd2010dec, cwd2010feb, cwd2010jan, cwd2010jul, cwd2010jun, cwd2010mar, cwd2010may, cwd2010nov, cwd2010oct, cwd2010sep
## min values : 0, 0, 0, 0, 0, 0, 0, ?, ?, ?, ?, ?
## max values : 150.96939, 205.35623, 47.32375, 65.93313, 32.79562, 233.96126, 219.87500, ?, ?, ?, ?, ?
This allows us to do operations much easier. For instance we can just run arithmetic functions on the entire stack:
CWD_mn_2010 = mean(cwd_stack)
plot(CWD_mn_2010)
We can also do ‘global’ statistics like finding the ‘mean’ of all cells for any given month.
cellStats(cwd_stack, stat='mean', na.rm=TRUE)
## cwd2010apr cwd2010aug cwd2010dec cwd2010feb cwd2010jan cwd2010jul
## 45.5497388 157.9214946 1.6478612 5.2903873 0.1809682 176.2727036
## cwd2010jun cwd2010mar cwd2010may cwd2010nov cwd2010oct cwd2010sep
## 143.6983408 27.6639933 90.9696350 19.2020629 27.7737959 120.8165948
Or other summary visualization like histograms on a month by month basis.
hist(cwd_stack)
Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case raster cropping and masking are useful for unifying the spatial extent of input data. Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.
We will use two objects to illustrate raster cropping:
srtm representing elevation (meters
above sea level) in south-western Utahzion representing Zion National
ParkBoth target and cropping objects must have the same projection. The
following code chunk therefore not only reads the datasets from the
spDataLarge package. This package needs to be installed
directly from github using the following code:
install.packages("spDataLarge", repos = "https://geocompr.r-universe.dev")
We will also need to import the sf package to handle
vector data.
if(!require("spDataLarge")) install.packages("spDataLarge")
## Loading required package: spDataLarge
if(!require("sf")) install.packages("sf")
srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))
Let’s take a look at the two files, note the use of
add=TRUE in the second plot command. This adds the new
object to the previous plot.
plot(srtm)
plot(zion, add=TRUE)
We use crop() from the raster package to crop the srtm raster. The function reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument.
srtm_cropped = crop(srtm, zion)
plot(srtm_cropped)
plot(zion, add=TRUE)
Note that the extent of the
strm raster has been cropped to
fit the zion polygon.
Related to crop() is the terra function mask(), which sets values outside of the bounds of the object passed to its second argument to NA. The following command therefore masks every cell outside of the Zion National Park boundaries:
srtm_masked = mask(srtm, zion)
plot(srtm_masked)
Importantly, we want to use both crop() and mask() together in most cases. This combination of functions would (a) limit the raster’s extent to our area of interest and then (b) replace all of the values outside of the area to NA.
srtm_cropped = crop(srtm, zion)
srtm_final = mask(srtm_cropped, zion)
plot(srtm_final)
Changing the settings of mask() yields different results. Setting updatevalue = 0, for example, will set all pixels outside the national park to 0. Setting inverse = TRUE will mask everything inside the bounds of the park (see ?mask for details)
srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
plot(srtm_inv_masked)
Raster extraction is the process of identifying and returning the values associated with a ‘target’ raster at specific locations, based on a (typically vector) geographic ‘selector’ object. The results depend on the type of selector used (points, lines or polygons) and arguments passed to the extract() function, which we use to demonstrate raster extraction. The reverse of raster extraction — assigning raster cell values based on vector objects — is rasterization.
The basic example is of extracting the value of a raster cell at
specific points. For this purpose, we will use zion_points,
which contain a sample of 30 locations within the Zion National
Park.
data("zion_points", package = "spDataLarge")
plot(srtm)
plot(zion, add=TRUE, col=NA)
plot(zion_points, add=TRUE, col='red')
The following command extracts elevation values from
srtm and creates a data frame with points IDs
(one value per vector’s row) and related srtm values for
each point. Now, we can add the resulting object to our
zion_points dataset with the cbind()
function:
elevation = extract(srtm, zion_points)
zion_points = cbind(zion_points, elevation) # insert the elevation data into the orginal shapefile
zion_points
## Simple feature collection with 30 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2077 ymin: 37.16632 xmax: -112.8717 ymax: 37.43165
## Geodetic CRS: WGS 84
## First 10 features:
## elevation geometry
## 1 1802 POINT (-112.9159 37.20013)
## 2 2433 POINT (-113.0937 37.39263)
## 3 1886 POINT (-113.0246 37.33466)
## 4 1370 POINT (-112.9611 37.24326)
## 5 1452 POINT (-112.9898 37.20847)
## 6 1635 POINT (-112.8807 37.19319)
## 7 1380 POINT (-113.0505 37.24061)
## 8 2032 POINT (-113.0953 37.34965)
## 9 1830 POINT (-113.0362 37.31429)
## 10 1860 POINT (-113.2077 37.43165)
Note that now we have the elevation values for each one of the points
associated with zion_points. For fun let’s use
mapview to make an interactive map of both. Again we will
use zcol to set the column used for symbology in the point
map.
if(!require("mapview")) install.packages("mapview")
## Loading required package: mapview
mapview(srtm)+
mapview(zion_points, zcol='elevation')
Let’s tie this lesson back to our work with OSM. Here we are going to
import the package osmdata to allow us to interact with the
Overpass Turbo API. Let’s try to build a query using a bounding box from
getbb(), where we search for Zion National Park, this
bounding box then get passed to Overpass Turbo and restaurants will be
added.
if(!require("osmdata")) install.packages("osmdata")
query = getbb(place_name="Zion National Park, Utah, USA") %>%
opq() %>%
add_osm_feature(key="amenity", value="restaurant")
restaurant = osmdata_sf(q=query)
Ok, bad request! Thats not too specific. Let’s try to see what happened. First let’s see if we were able to get a bounding box for Zion National Park.
getbb(place_name="Zion National Park, Utah, USA")
Looks like that is the problem. You can try to see if you can adjust
the name to get a valid response. But its likely we will need to create
a bounding box manually. Luckily we can get the bounding box from the
srtm raster file using bbox(srtm), this then
needs to be converted to a vector using c() in order to
meet the requirement of the opq function. Look at
?opq, notice that we can pass a bounding box in the format
bbox = c(xmin, ymin, xmax, ymax).
Note: bounding boxes can also be found using tools like bboxfinder.com.
Now let’s try again, using our new bounding box.
query = opq(bbox = c(bbox(srtm))) %>%
add_osm_feature(key="amenity", value="restaurant")
restaurant = osmdata_sf(q=query)
restaurant
## Object of class 'osmdata' with:
## $bbox : 37.1320834298579,-113.239583212784,37.5129167631658,-112.85208321281
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 32 points
## $osm_lines : NULL
## $osm_polygons : 'sf' Simple Features Collection with 4 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULL
Now let’s plot it out and assign the restaurant names as labels:
mapview(srtm)+
mapview(restaurant$osm_points,
label = restaurant$osm_points$name,
legend=F)
Raster extraction also works with line selectors. Then, it extracts one value for each raster cell touched by a line. In this case, the best approach is to split the line into many points and then extract the values for these points. To demonstrate this, the code below creates zion_transect, a straight line going from northwest to southeast of the Zion National Park
Let’s create a line based on two points
c(-113.2, -112.9) and c(37.45, 37.2) which are
xy pairs in lat lon. Note that we are assigning a projection using an
EPSG code of 4326. EPSG codes are unique idenitifying codes
that are availabel for all common projections. To look up EPSG codes
please go to https://epsg.io.
zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
st_linestring() %>%
st_sfc(crs = 'EPSG:4326')%>%
st_sf()
zion_transect
## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## geometry
## 1 LINESTRING (-113.2 37.45, -...
Now let’s see if its showing up in the right place
mapview(srtm)+
mapview(zion_transect, color='white')
The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an ‘elevation profile’ of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.
Because our line is made up of only two points and
extract works best with points, we need to break our longer
line into lots of points.
The first step is to add a unique id for each transect. Next, with
the st_segmentize() function we can add points along our
line(s) with a provided density (dfMaxLength) and convert
them into points with st_cast().
Let’s look at st_segmentize help to see what units
dfMaxLength is in:
?st_segmentize
Now lets create our collection of points:
zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")
zion_transect
## Simple feature collection with 257 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## First 10 features:
## id geometry
## 1 1 POINT (-113.2 37.45)
## 1.1 1 POINT (-113.1988 37.44902)
## 1.2 1 POINT (-113.1976 37.44805)
## 1.3 1 POINT (-113.1965 37.44707)
## 1.4 1 POINT (-113.1953 37.4461)
## 1.5 1 POINT (-113.1941 37.44512)
## 1.6 1 POINT (-113.1929 37.44415)
## 1.7 1 POINT (-113.1918 37.44317)
## 1.8 1 POINT (-113.1906 37.4422)
## 1.9 1 POINT (-113.1894 37.44122)
Now, we have a large set of points, and we want to derive a distance between the first point in our transects and each of the subsequent points. In this case, we only have one transect, but the code, in principle, should work on any number of transects:
zion_transect = zion_transect %>%
group_by(id) %>%
mutate(dist = st_distance(geometry)[, 1])
zion_transect
## Simple feature collection with 257 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## # A tibble: 257 × 3
## # Groups: id [1]
## id geometry dist
## * <int> <POINT [°]> [m]
## 1 1 (-113.2 37.45) 0
## 2 1 (-113.1988 37.44902) 150.
## 3 1 (-113.1976 37.44805) 300.
## 4 1 (-113.1965 37.44707) 450.
## 5 1 (-113.1953 37.4461) 600.
## 6 1 (-113.1941 37.44512) 750.
## 7 1 (-113.1929 37.44415) 901.
## 8 1 (-113.1918 37.44317) 1051.
## 9 1 (-113.1906 37.4422) 1201.
## 10 1 (-113.1894 37.44122) 1351.
## # … with 247 more rows
Finally, we can extract elevation values for each point in our transects and combine this information with our main object.
zion_elev = extract(srtm, zion_transect)
zion_transect = cbind(zion_transect, zion_elev)
zion_transect
## Simple feature collection with 257 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS: WGS 84
## First 10 features:
## id dist zion_elev geometry
## 1 1 0.0000 [m] 2001 POINT (-113.2 37.45)
## 2 1 150.0961 [m] 2033 POINT (-113.1988 37.44902)
## 3 1 300.1922 [m] 2020 POINT (-113.1976 37.44805)
## 4 1 450.2883 [m] 1989 POINT (-113.1965 37.44707)
## 5 1 600.3844 [m] 1902 POINT (-113.1953 37.4461)
## 6 1 750.4805 [m] 1869 POINT (-113.1941 37.44512)
## 7 1 900.5766 [m] 1827 POINT (-113.1929 37.44415)
## 8 1 1050.6728 [m] 1767 POINT (-113.1918 37.44317)
## 9 1 1200.7689 [m] 1772 POINT (-113.1906 37.4422)
## 10 1 1350.8650 [m] 1722 POINT (-113.1894 37.44122)
We can now see that we have both the distance from the first point in our line, and the elevation of each point along that line.
Let’s plot the transect and see what it looks like:
plot(x=zion_transect$dist,
y = zion_transect$zion_elev,
type='l',
col='darkred',
main='Elevation Profile',
xlab = 'Distance',
ylab = 'Elevation (m)')
Step 1 - Create a query to pull all the ‘primary’ roads in Zion
National Park and save them in a variable called roads.
Make sure to look at the tagging convention on OSM’s map
features directory. Also you need to set the bbox from
the query, you will want to do it based on the srtm raster
as we did earlier.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# # hint
# query = opq(bbox = ?????) %>%
# add_osm_feature(key="?????", value="?????")
#
# roads = osmdata_sf(q=query)
Step 2: osmdata_sf returns points
(roads$osm_points), lines (roads$osm_lines),
and polygons(roads$osm_polygons). Create a new variable
called primary_roads that holds the road line features.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# print(primary_roads)
Step 3: Use the filter() function from
dpylr to isolate only roads with the name
Zion Park Boulevard, and store the output in a variable
called park_road:
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# print(park_road)
Let’s take a look at the data:
park_road
Notice that there are multiple segments making up this one road. Roads are typically broken into segments to allow for connectivity between other roads, and to maintain data like speed limits - which likely vary along the road.
Step 4: We need to ‘dissolve’ all of these road features into a
single line representing ‘Zion Park Boulevard’. In order to do this we
are going to tell the computer to create a new geometry
column using summarize, based on the st_union
of the multiple geometry rows. This is one that we haven’t talked about
yet, so I will just do it for you.
park_road = park_road %>%
summarize(geometry = st_union(geometry))
plot(park_road)
print(park_road)
## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -113.0368 ymin: 37.16021 xmax: -112.9902 ymax: 37.19982
## Geodetic CRS: WGS 84
## geometry
## 1 LINESTRING (-113.0368 37.16...
Notice there is now only one row of data and our multiple road segments are now the union.
Step 5: Now we can create our column called id, then use
st_segmentize to break the road into multiple points at a
fixed distance (250 or 500 works well), then we need to cast that point
object back into an sf point object using
st_cast.
See if you can find the correct chunk of code
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# park_road$id = 1:nrow(???)
# park_road = st_segmentize(park_road, dfMaxLength = ????)
# park_road = st_cast(park_road, "?????")
# park_road
Step 6: Let’s calculate the distances from the first point to all following points:
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# park_road = park_road %>%
# group_by(?????) %>%
# mutate(dist = st_distance(?????)[, 1])
# park_road
Step 7: Now we can use the extract function to grab the
values from the raster srtm at each one of the points
making up park_road, then we need to cbind
those values back into the park_road_elev dataset.
### ### ### ### ###
# YOUR CODE HERE
### ### ### ### ###
# park_road_elev = extract(????, park_road)
# park_road_elev = cbind(????, ????)
# park_road_elev
Step 8: Plot the park_road_elev$dist on the x axis and
park_road_elev$park_road_elev on the y axis, and set the
type to 'line'.